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Summary 

The Weibull distribution has been widely adopted for the statistical description and inference of 
fatigue data. This document provides user instructions, examples, and verification for software to analyze 
gear fatigue test data. The software was developed presuming the data are adequately modeled using a 
two-parameter Weibull distribution. The calculations arc based on likelihood methods, and the approach 
taken is valid for data that include type 1 censoring. The software was verified by reproducing results pub- 
lished by others. 


Introduction 

This document provides user instructions, examples, and verification data for software to analyze gear 
fatigue teat data. The software was developed on the presumption that the data are adequately modeled 
using a two-parameter Weihull distribution, it is based on the likelihood methods described by Meeker 
and Escobar ( 1998). The software can be used to determine 

I Maximum likelihood estimates of the Weibull distribution 

2. Data for contour plots of relative likelihood for two parameters 

3. Data for contour plots of joint confidence regions 

4. Data for the profile likelihood of the Weibull distribution parameters 

5. Data for the profile likelihood of any percentile of the distribution 

6. Likelihood- based confidence interv als for parameters and'or percentiles of the distribution 

The software can account for tests that are suspended without failure. The statistical terminology 
for suspended tests is "censoring.” The analytical approach for the softw are is valid for type 1 censoring, 
which is the removal of unfaikd units at a prespecified lime. Confidence regions and intervals are calcu- 
lated using the likelihood ratio method. Guidance for the philosophy and interpretation of statistical confi- 
dence interv als can be found in the text of Hahn and Meeker (1991 1 . 

The Weibull distribution has been widely adopted for the statistical description and inference of 
fatigue data. The Weibull cumulative distribution function F is defined as 


F< 


M* I -expi 



where t is the tunc lo failure. 0 is the shape parameter, and n is the scale parameter For the Weibull 
distribution, the probability density function can have u variety of shapes, and the hazard function may- 
be an increasing or decreasing function of time. The shapes of these curves depend only on the shape 
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parameter. Typically, surface fatigue data for gear* have distributions with shape parameter* less than 
three. Therefore, the probability distribution functions typically arc skewed right. For a given Weibull 
distribution, one can determine the mode, mean, median, variance, and higher moments analytically by 
employing the gamma function (Cohen. 1973 >. 

The software developed herein was written in the FORTRAN 90 programming language. Users of the 
software ate required to write short main programs that provide the data to be analyzed and also call the 
appropriate subroutines. The examples provided in this manual were compiled using the Microsoft 
PowcrStation Fortran Compiler, version 4,0 (Professional F.dition). 

The data analyzed as examples for this document (tables 1 and II) are those from gear fatigue experi- 
ments. These data have previously been analyzed using regression-bused methods, and the results 
reported ( Kmttz et al.. 2000). The data are again analyzed and reported in this document using software 
based on likelihood methods. 

Verification information for the software is provided in appendix A. The softw arc was verified by 
reproducing results published by Meeker and F.scohar (19V#). Appendix B prov ides descriptions and 
instructions for using subroutines. Some calculations require the use of routines from a standard math 
subroutine library' flMSL. 1997; IMSL is a registered trademark of Visual Numerics. Inc.) or an equiva- 
lent Subroutine »maxU (appendix B) i* based on the code of Keats. Lawrence, and Wang ( 1997 >. 


Estimates of Distribution Parameters 

DataSet 1 

The main program that analyzed the data of table I to provide maximum likelihood estimates of the 
Weibull distribution shape and scale parameters is presented below a* source code listing I Allocatable 
dimensioning is used (lines 8 and 9* to make certain that the array size* match the value of “n. ’ the vari- 
able describing the number of lests (line 7). For the army • censor.'’ a value of "zero" is given to indicate a 
test run to failure whereas a value of "one'' is given to indicate a test suspended without failure (lines 1 0 
10 15i. A call to subroutine h-rrexr// provides maximum likelihood estimates of shape and scale parameter* 
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(tine 40) With (he maximum likelihood vahte* for the shape and scale available, subroutine wblm-erre 
provide* maximum likelihood estimates for any desired percentile of interest (lines 48 to 52). 

The main prof rim produced the output presented below following source code listing 1. The maxi- 
mum likelihood estimates for shape and scale parameters are 1 .037 and 376.7, respectively. 


IW f*r» tr.+ lir.e Kiciwtoft Ter v nr. Cvn >****• 4.0 


2 . r*if. pidgrat 9 * »r *.«r. rtfiUi 

? tip lie it t*/r.* 

> irtr^pr r . . i 
4 r**l. «1 *» 

« irrrflpr, *11 rvat.ahie t • ccn»st':1 
4 ml «***.«, *r»r . 1 if* 

t n n 14 

* #V <r.t ) 

* alloct** ' w,*cr p. ' 
iC dc t * . , t 

ii ctnocrl * 5 
*2 •r.f ds 
W 4c ;.9.r» 

14 c«n»cr*k: • 1 
14 *r.2 ds 

14 cprt' 2 ; (u»**i 4 «rlK 2 txt* 

;*> 1 mi rt mtt tc tms ot tsst Tints 

:« UM 1 • 45.51 

11 tin* 5 * 51.44 

55 UM 3 ■ 54.14 

s; um 4 • <mi 

52 UM 5 • 15.74 

51 tint < • 1 14 52 

54 tit* 7 . 1H 41 

5* tlM I • 254.25 

24 tin# 9 • 299 55 

57 t;*c iC - *55 54 

54 t^M 'll - 521 44 

29 tlV I? • 124 21 

>2 n • 1 7* 

i: tit* 14 1* )V K 

32 t'M* :5 • >15 4f 

ii 4c i*;.n 

44 if rc#r.«cr'l< c > t**r 

i< «1 m if ictraord' *e .i t.vfl 
3“> writ* '20 . 391 * 

21 end if 
29 «r>d do 

45 c*2i wr**!; tit* ronoor.n. 

42 fcrit* *29,225' 

42 *rU*'20.i;>C' *t*?* 

41 520 for~*t I IC .5 . * test iB ii«d •) 

44 ill Ion«t f 1C .5 . * t*4t 

44 125 fomt ’ l*fc«5lhco4 MUMt*i * 

44 : >4 (orMt ‘ • ',*12.4 ‘ »:kit • ‘.4l2 4 • 

4^ frit • 2 - C 
4 1 do t • 1 . * 

44 fr*c • fra* • ft.i 

5 *j c*i‘. *fcilr./«x »• tBn*?* IfOt . l-f« 

11 *ri?#'25. ii? fr*c *ifa 
M tn4 dr 

41 .77 f« 1 fmj '.rt -or. Ilf* • * ,flC 4 

44 rlCp«r'25 

55 **05 
M «:4 
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Output produced by this main program waa a» follows: 


4C.2( ;or. 

n.M ttit 

on t««t f#i;#* 
imi t»u 

0^1 ttit Inin 
1);.22 test t#i ltd 
omi isit u.uj 
; r >' 2* test f#.ltd 
?l* ^ test luflpsrxliC 
K2.*f tect iueptndtd 
10; ** T*ftf i^tpended 
104 0» t#*r tu*p*n<4«4 
20* • ! test •utptrvJed 
20« S», test •u*p«n4*5 
test •utptndtd 

mmxr jm lifctliboo* estlnetes 
shefe * ecele • .H57I03 

. ZC pxc portion lift • 4;.C4CJ» 

■ 2( proportion lift * •(.‘’ill 

K proportion life • IH.C* 05 
*C propoTtirn Irf* * 1*1.1*52 

1C prepext icn life • 2*4 €C' 4 

(C prcpcxrien life * U€ r 
. **£ proportion life ■ *5C 5C«9 

. 1C pxspoxticn life * 5K.MC1 
. >C proportion life • til. (HI 


Esimatti of Distribution Parameter* 

Data Set 2 


The main program that analyzed the data of table II to provide maximum likelihood estimate* of the 
Weibull dictrihutioo shape and scale parameters i* presented below a* source code listing 2. The source 
code listing 2 minors source code listing I. differing only in the lines defining the number and results of 
experiment* (lines 7 to 32). The main program produced the output presented below' following soutce 
code listing 2. The maximum likelihood estimates for shape and scale parameters are 0.8616 and 103.0. 
respectively. 



i . as-.r. pro*™.- te anal »jp*r fl***f u*»r ter-.^r ♦*«» result* 
1 ^litH mm 
1 .nlt^tr u.m 

i it#.. tiltottt&it. t ttr#' 

5 , g . locat ao 1§ *. i ctr.#07(.i 

4 r««l «Mpt. tx*;,litt 
- n * 2C 

k il’sc^tc 'time a ; 

* ftKec«r*'ctn#:;r t a) • 

do i«i,n 

• c 

;; rr< <%r 
\t 6c i*:t.r 
24 x#nt©r;i • : 

1! er 4 ffc 

1< optr >: fil#»-#Jp#rrt#?l txt # > 

n . he* i nr **r% ro »m* the a**ay or urr t:*u 

;» vc« 

it I n 
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ac 

tt«t 

2 • 15.12 


a* 

Uftf 

4 • 3 2.24 


22 

ti«C 

> • JI.2'' 


;j 

ttu 

4 • 25. n 


at 

ttu 

r- 21.22 


as 

C.M 

t • 40.21 


at 

t;w 

*:• 40. *>* 


a* 

tit« 

!$>• 41.42 


at 

t ;«• 

i;>* si ?r. 


as 

t l*c 

»a2»» 5*. as 


v» 

f l*** 

<)!;• 54 ?*> 


i; 

t i%* 

• :* • ?c . *: 


>2 

t in* 

-IS*. 12*>.C'’ 


13 

dir* 

«14i* 132.40 


24 

Sir* 

‘i’i* 133.31 


i/ 

Or# 

--0.02 


>4 

Or-* 

•2»Jy aCd.Od 


a" 

tn 

(22». 2CD.02 



*0 S' 

• ; ,r. 


24 

a 

»c*:>*or{ii .*q. 

C ; t.*n 

if 

rrit< 

4 2C.1CQI tlM'i: 


4 C 

• 1** 

if <ct^*or(.i . 

•q. 1 ) than 

a 

»;.v 

• 1C. 1C*.) t iff.* ' 1 1 


«a 

tr.l 

*f 


42 

•r«* 

tr> 


44 

call 

inuxli t i^e , senior , n. *hap* , e 

45 

writ 

# ac.ics. 


44 

writ 

e 20,104* aMr*. 

•nlr 

4? 

ICO 

for-at (f 1 0. 2, * teat 4ttl#4 ') 

41 

1C) 

for-*t if 1 0 .2 , f ttrtt mufervSc* *) 

4 9 

105 

format ** Mxmr 

1 t/r) ihcod eaMfta 1 

SC 

io* 

format ' * afwipe • 

* ( r:2.4 , * ucal* * 

%) 

!rr 

• o.r 


52 

6c i 

*i » 


n 

for 

• f r« r * C ; 


*4 

call 

wklir.var** raha?t . act** . f rac.lift. 

•5 

writ# 2 d, idd, fr*c lift 

SC 

#rv5 

6 c. 



5 ^ id' fcr**t ft.;,' 4 report icc lift • *,t;:.4 
!l clc.at'2d» 

£9 ttCf 
42 r>i 


Output produced by this main program was as follows: 

''.Cl cut (uUd 
ctit fiii*4 
ii Oi tut £*il*i 

a; 24 tut 

>3.?*? c*i*. fail** 

is ic t*at fail** 

31 j: ttit U;U5 
4c:s t*«t 
ic . * 4 teat fail** 
iM*» tilt uUtj 
Vi. '*C t*at (i;Ut 
jt ;> tut u.ui 
it.i 1 tut fitli4 
:c.i; tut UUil 
is’j.c'’ tut tail** 

1IC.S5 tilt tllilj 
;n.3i tut u;;h 
KC.C: tut awipart** 
iCC.05 tut » 

ICC. 0 V tut *4 •**.10*4 

;ihtUhXi llti'ltti 

IM^ • .Itut«cc teal* • ,:-4l0t*0*» 

.Id pr&^ttur. lift • JJ ? i 
,2i (. report k tr. iif* • 1I.JS54 

.1! f report ic,r. lifi • 11 Hid 
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.«C fropcrtlcn 
. tc fr?if«rt*sn 
proper* Jr.-. 
. *c proper* Jen 
•0 pro fvr* ion 
. fC prcp&rti»rv 


uu 

0 

r.aijs 

lift 

0 


lift 

0 


;jft 

0 

tsr *;cic 

: \ 

0 




2 tm .<445 


Joint Confidence Regions 

DataSet 1 

The main program that analyzed the data of table I to determine joint confidence regions for the 
simultaneous estimates of the shape parameter and 10-percent life is presented below as source code list- 
ing 3. Output files were written to provide data for likelihood ratios and confidence levels (lines 25 to 28). 
Allocatabie dimensioning is used (lines 30 and 3 1 ) to make certain that the array sizes match the value of 
"n," the variable describing the number of tests (line 29). For the array ‘'censor.” a value of “zero” is 
given to indicate a test run to failure whereas a value of “one” is given to indicate a test suspended with- 
out failure (lines 32 to 37). The range of values for the shape and scale parameters must be provided (lines 
59 to 62). A cal) to subroutine \ jlotue provides the data required for contour plots. The subroutine writes 
results to FORTRAN' units numbered 9 and 10. The output from the program is two data files, each con- 
taining 1681 lines. Each line consists of three numbers corresponding to a 10-percent life value, a shape 
parameter value, and either the likelihood ratio or the likelihood-ratio-based joint confidence value. The 
joint confidence values are displayed as a contour plot (fig. I ). providing a graphical display of plausible 
values for simultaneous estimates of the 10-percent life and the shape parameter. 

7 
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hr* fcoerre Li re * 1 — naof t fcrttar Pouarfitat m. C©-pi-er. Version 4.0 

; . * Min prc^re^ to calculate likelihood ba**d prcbeM **.*y cent our a 
; and relative lihei-aocd c inters ©! WeihuH distritut ;or.s fro* censored date 

i ; 

4 . tni* pro gra« p~f*r<u-*a that' the two pax amt tars t© 

i ; fc* plot?** are the shape paramatar and '.ta 251 lift 
f • version 10 

1 written fcy Ts« Kran*f Or. 2*12-2052 
f impliM* non# . declare 412 variables 

5 Intent r n . the number ci testa 

;c raal . allocetatle : : tie*:: : the tart times 

21 integer. ailocstafcie : eertor(< ; th* ceraonng inf ©mat. or 

i; raal terra* . tetwin s the max and m*n vnlu#* for th? 15 par-cant Ufa 
;i » to te ua#d in calculation* 

24 raal shaper in shapemax tha ra* and min valve* for the ahapn f art or 
15 i to ta used in cal cular ions 

14 integer i . a do loop countar 

15 i am aa needed for a particular cat# 

If 1 c;**nga only line* ©etween tha two lines cf atari 

1» l 

2t i tV outputs for th:* prcgr*~ w-22 be written tc thtaa two frit* 

22 : the file format?. art a© follows <aach he* 1 coljmna 

22 : coVimr 2 ia to# shape parameter 

21 . coiumr ; ia tha 1C pc ream Ufa 

24 . fclumr 1 la tv li>.elihesd rati© cr tha serf ider.ee numfcer aa a percentage 

it . file for -jr.tf #:f ia rhe confidence r.jtlic data 

24 : flit for unit 19 :% the likelihood ratio data 

29 0Pfa: -C. IILF * *euper proh.tr") 

2i OffiV > KiL» • •rep'*’- • rar \n t x?'t 

if n • 23 

iC allocate time 'ni . 

>; a:;oc«r.« ceneor tn) • 

>2 dc 1*2.1 
11 carasr'i; * c 

14 and do 

15 do 1 • I . r. 

14 canacrU) » 2 
19 #r^ do 

it : xs*t k* inris tc tyrir thk a*j*y cr rn: ?:rcs 

it tie* 1 • 43.24CC 

40 UM 2 * U.MCO 
42 Kiw* 1 * M MC9 
42 tin* 4 - fVMCS 
42 tir* 5 • >i 7 %t0 

44 tiw < • 114.223 

45 tif4-5 • lit. 413 

44 Ufa I - 254 2C0 
49 tire *1 » 21* 320 

41 V.at 1C - 221 540 

45 time 11 • 251.440 
EC U«* 12 • 234 0*0 

51 h— II - 255 ICC 

52 t.me 14 • 254 ICO 
f) t.me 15 • 2I5.4CC 

54 • here »e #at 1 1* t nt.U ioz 1C ptrvtr.t liil and ahapa factors 

55 * tenrax as th« lasgaat 15 parcant lifa valua 
?< : tfiwin is tha sataUae: 13 partar.t 1-fa valua 
a 1 iiMp#*Af is the largaa? lU^t val-t 

•I ah^'-1r ii tha s-a licet shape value 
•4 tenray ■ i«' 

4 ; t#rw;r. • i t 
4. anafm^ax * 2 

42 shafamir • *. 4 

l j r ^a sf tha lines be '. 0 * na#d fn h* rrarged tc run a part ic\uar 

45 the call tc this subrout ire for arrort i r. the rancor iryj srxa/ 

44 call cchec* fr.. censor » 

4 7 . new wa asit for tha solution. 

41 call time, censor, r, shape-ax ahener-.n ter - ax, *#nnlr.» 
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The main program (Km analyzed the dau of table II to determine joint confidence region* for the 
simultaneous cm i mates of the shape parameter and JO- percent life is presented below as source code list- 
ing 4. The source code listing 4 minors source code listing 3. The output from the program is two daw 
files, each containing 1681 lines* Each line consists of three numbers corresponding to a 10-perccnt life 
value, a shape parameter value, and either the likelihood ratio or the likelihood-ratio-based joint confi- 
dence value* The joint confidence values are displayed as a contour plot (fig. 2), providing plausible 
values for simultaneous estimates of the 10-percent life and the shape parameter. 



Li »t »n2.* 


Lin* Scarce iin* Macic# s»fi f* rti'Sfi Cs^ilir Vr,,rf 4 


2 

1 

4 

I 


Vi *•••* *rc<#bi-lry rri-cjrr 

rf! rrUtlv# UK.lihooJ vv.V.c-i* »t Utsfeui; 4, «t .•»**»<*• tr~ ««errrd <Ut. 


th;.* sroyr*# p;titrii t;«*t tl*« t^c lo 

b* ir« to* tA*s« to# I - "- ** ft 
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4 r* rater, i c 

? aritttn Iff 7 ,m tr*r*t t, f 2 : W'.C; 

4 implicit non# r *#r:ar* all viriillti 

I inttoir n tht mithtr r > f tent* 

it rtAi, Alloeatafci# 5 . the t«at "i*ea 

ii ir.trjtr •::oc»t#h> : etr.nor « - * r.h* ctnocrir.? inf or rat .or. 

.2 rtAl ttr.***, ttn«ir> t th* max av! tin values for rht 19 p*rctf.. iift 
ij te b« ward in calculation# 

14 rtal BhApruin »h«perAx 1 tV «uix ani -ir. v&Vje.e for the chape factor 
;« • to be u»tA In raieuUriors 

it ir.tactr ; i a do loop ecur.ttr 

1? . cnarvjt proqra* ma nt*<Vi for « part* ot*. ar salt* 

11 . entry* only lint* bttwten tv two lmtn of §tars 

if ; * ••••* 

2C I tht output* for this p?5<?r*» will bo writtfr tr these rwo filer 
;; 1 tht flit format* are ai f el lew* :e*eh Mr 1 ml aatsi 

12 t sol'jwr. 11* tht «h«pt parameter 

21 l coljnn 2 1 * tht 10 percent lift 

24 I coitw 2 U ttt UfcAilhoo* r*t ie or t be cstMdfVf r.ut*.*r ** n peroerttgr 

25 . file for jr.it lie it the confidence mmhtr da- a 

24 ; fiit for jr.it • * l* tht Uktiihcod ratio data 

2*) OF CM 19, Fil.* - *fca*t*profc txt*» 

it GKK I. Pii4 * •bt*t*rttvo txv * 

24 fi • 22 

1C ailocatt ’t-e* r. j 
K ailocata ct.tatr ra- 


12 

Ac i* 

1 . I** 



n 

ceaai; 

r 'i » 

• C 


14 

t»A >o 



IS 

Ac 1* 

21. r. 



14 

carta i 

r 'i 

• • 


1- 

tr^A Ac 



11 

( Kfft? « 

tme to 

If 

UM 

2 ) » 

•\o« 

49 

UM 

2>* 


S 

4. 

titMr 

ft>« 

is 

52 

42 

t.rf 

' 4 ) ■ 

21 

14 

42 

tire 

fftj* 

13. 

1- 

44 

tic*# 

•4; • 

1! 


45 

t*C4 


21 

22 

44 

t ir<t 

'•»* 

4 1 

24 

4** 

txre 


4C. 

**4 

41 

t i*t 

10» • 

4ft. 

4? 

4» 

lift* 

12> • 

s: 

*C 

SC 

tia* 

12)* 

55 

2> 

s; 

tin 

n< 

5t 

A f 

S3 

ti*t 

14)* 

?t. 

i; 

S3 

tm 

IS- • 

12 1 

.C“ 

14 

tint 

11. * 

123 

.15 

!? 

tint 

n>. 

111 

.21 

54 

ti*t 

11;* 

JCC 

i.CC 

ft ’ 

tl*t 

l*t* 

j33 

C3 

SI 

tier 

2 3) • 

>05 

i.CC 


** . htr* vt Ht tht Pernod* frr rerrar.t nr# ana »nap" 

45 te?*a* i* tht Ur9*t* 1' P*rrtn* !!•* vatur 
4; tone in 1* tht rrallt*c 10 percent lif* value 
4^ *hapt**x tt tht larqta* *h*f* valu^ 

4 > *hapr**.n i* tv mlUv ah ape val i» 

a « ttr-ax • ) ;> 

49 ter-ir. . ;.i 
44 i:vyri» • . 4 
4 41 an*p «*:r. * 0 4 

44 rent of tht ;ir*a Mi*.* n#*S to ht efcar.s** to rjr. a psrt.euUi caac 


V, tht call to tM* a ihrs *Mn# *f.r trrorr ir "he ctranrirg array 

oar. erherV*n.e#f.ehr; 

*»2 r.rv -at a** for tv *c'j?‘.or 

s, eai: w*lcfvt:ti*t rtraor.r. Aft*pa**v aiv.pe«in rer-ax, *asitin 
n ?ic*t 
*»i eloet'f* 
u atop 

n- tfji 
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Profile L&eiihood and Proflie»I JkeUhood-Based Confidence Interval* 

DataSet 1 


The main program that determined the profile likelihood and profilc-likelibood-based confidence 
interval* for the 10-percent life for the data of tabic I is presented below as source code listing 5. The 
mam program determines the profile likelihood over a range of value* for a specified percentile of inter- 
est. m this case the 10th percentile. Any percentile can be investigated by changing one line in the code 
for variable "pfrac" (line 49). Aliocotablc dimensioning is used (lines lb and 19) to make certain that the 
array sires match the value of "n." the variable describing the number of tests (line 1 7 ). For the amty 
"censor.’' a value of "wo" is given to indicate a test run to failure whereas a value of "one" is given to 
indicate a test suspended without failure (lines 24 to 29). Array "time" (lines 31 to 45) provides the test 
times. The subroutine ptvpercent is called (line 56) to determine the relative likelihood ratio for a speci- 
fied 10-percent life. By calling the subroutine repeatedly within a do-loop (lines 50 to 59). data for a 
profile-likelihood plot art calculated. The output from the main program was plotted using two vertical 
scales to show both the profile-likelihood ratio values and the confidence level (fig. 3). Table ID lists 
some profile-likelihood ratios and the corresponding confidence level. A 90-percent confidence interval 
for the 10-petcent life is shown graphically on figure 3. The low er and upper bounds can be determined 
more precisely by running the program again with an appropriate range of values for the 10- percent life 
(figs 4 and 5). 

tab: s »i —confidence levels 

AMI fOMtf.SaerNDING PROFILE- 


I IKK IHOOO RATIOS 



10-par ca ot Ma. mMona of cyctaa 

Figure */— #ro(Ua-tfco*bood raooa and conManca tavat# tor 10-pafeant Ma for data of 
rental. 
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1»-parc*nt Ufa, mfMone of eyetot 

^ <■— Pro#t«-W^nood-ba**d conftdanca Interval for 10-percent lit* tor data of 
tab)* I (tow or bound). 


0.30 — 



npuo 9,— fmdto-Kfceltoood'toaeed confidence totorv*) tor lO-percent If* lor data of 
tabfa I (upper bound). 


The ssibroutinc ph^ffnenr called by the main program determine* the shape and scale factors corre- 
sponding to the largest likelihood value by an iterative search method. This method requires that maxi- 
mum and minimum values for the shape parameter he provided to begin the search. The minimum and 
maximum values used in this example were hardcoded in the subroutine as 0.03 and 80 respectively. 
Users who with to analyze data sets that require inspection of shape vulucs outside the just-stated range 
need to modify the subroutine (appendix Bj. 


tur;» IZii iuiixq i 


-lr« tevret -;rw nctm ft ter :r*r. ►wtrfct At lor v«r»*9n i.t 

• i a •‘•if, £ x 07 r am i& c«.:.*At» Nh; 

2 . lakthx; j y«i of W**Lv«l 4ealr**>wl*&r*« figt. c;«uvr«'J j«ti 
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4 t pi *».-*• ft pftTtftf,*. 

; . i if* :.>♦ :f 

f . truiiar T.r K*a.v.r ftr. 2 ft a *2 ft. l 

* 4»{. lit ; *. r. j/r.a f it*:**rt •«» vftrifttl** 

> tn* (>:.9kr. rtprftttn*. ir.fJtft 
< . •““«*/ «. * Utd ;ow«4 ».t tin* »4.r* ( rvjr •* 
jr sr.*.«firi r .* t-'vt nu-ft^r c # ?<*4‘** 

*1 *1 , vnut!.* ! * tM tft**: ti**i 

;? ;p. «] ft 'i'-Ar : : rer.i^i : i t.«« ^antprin? inf rrmft*. i in 

i*it*3*r ‘ i dO .V.p CSJT.tt r 

•4 r«*l Sft* ftptl • 

;s -tftl pfrfto 
it rftft. ciran.t^ax 4\ 

!*» r. • 15 

:« *;io:ftt* o i • p. • 
ft* ft.iGTfttt center *r. > 

30 Ktft£ VZ VI tZ TC CkUJOKHC TKF^rATm 

21 A VAiLt C> 2* KiA?>S THAT A ffcTV«f 

22 . A VALVfc TF •-** THAT TEST VAT :TIi:SOfElft KITH VZ FAJL’JFE 

2> . KTI* 7>l VAAUtt A4 irTS0E?C 

24 i* ;*1 t 

25 • C 
2< ft-'ifc Oo 

2* do 

21 oer.tcr ' ; • 2 

ft* trA «k. 

v vrrF kk ntf- r: trrrf m AMA7 or tilt ;w 

3 ; 4 c 2<TC 

>; • *? %<ro 

j; *.* itr o 

>4 

it U*4 5>» :•! "COO 
*4 0 • 134.220 

i'' 5i*4 ' lit 4ft j 

it t^M ft** 25 4.502 
ii t ft * i ft >. *2: 

4C * *•*»•; • it.. 54; 

41 ♦.i'M .l » 20 C4> 

42 ti*M * .a * iC4 eft > 

4 2 t’."* ' *i • 20 . •'> ft 
44 t ;"4 .4 • 2C£ . 12 ft 
4 it ti** ft 5 * >3 0. 4 ft ft 

44 • ?h* c«:i tc thift eh^fffcft for tuvx* -n t.^ ••nftoria* ftrrfty 

4^ calft c?Y*rYi -r.rer.tftr 

41 op*r. it.FIif* ‘r^r^K J&rr.C lift .tut- 

4ft p fr«c * ft a 

•ft t^ir • 4 

*; t m *T. • ;30 

52 do • »?*t* r»*.p.i ' ftT,o. 

f> t*ir • r^ir d* 

54 6c ;*i . *01 
5ft t»ir. • twin ♦ 

54 call p: *p+r~+n* •r.m w.wrj.. ra^.fciv.apt •*/ tflfc./ 

5^ vrlt#a2»l r S' T»:n.«r hl v * D * v 

54 .5; for^Aft 4*ft0.T 
5ft trvd 4e 
I. t;op 
4; •ni 


Data Set 2 

The main program that determined the profile likelihood and profile-likelihood-baaed I confidence 
iMeivalt for the 10- percent life for the data of table II i» pte tented below at tource code listing 6. The 
main program mirror* that of tourcc code 5 The output from the main program wa* plotted using 
two vertical scales to show both the profile-likelihood ratio value* and the confidence level (fig. c). 

A 90-percent confidence interval for the 10-percem life u thown gmphicall> on figure 6. 
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o 2 « « • 10 12 14 1« 1« 20 

10-porcant ¥•. mMions 0 t cyctss 

#. Prott* NMtnood ratios and confldane* lavat* tor 10-parcant tt* tor data of 
tabtall. 


The subroutine phperctnt called by the main program determines the shape and scale factors corre- 
sponding to the largest likelihood value by an iterative search method. This method requires that maxi- 
mum and minimum values for the shape parameter be provided to begin the search. The minimum and 
maximum values used in this example were hardcoded in the subroutine as 0.05 and 8.0 respectively. 
Users who wish to analyze data sets that require inspection of shape values outside the just-stated range 
need to modify the subroutine (appendix B). 


SSxlSS ai&LUa s 

Scare* a:a* r.ttrtuft nrxr *r. Scv-rf- v v*r*,w. '..c 

1 • *r»ir. pre^r*.- to calejlsr* ' ISalir.oot 

2 ■ JaKtlioo^d curfti c t 41**?: rr»? 

\ . tizn c#r.*ort3 

t . 

5 ' tti* pro;r*r $>z U**t a p*r:«rt.;r if ty or of 

< • vrn l r.T, ; t 

^ ry T;f* r**r.;a On 2*13-233; 

I l-p* *. -if r,fv . v*r 

> ! th# f;11r.vlr.| v*ri*n >3 rvpriMAt infrr*it» 

10 * thty «r# JuM-crvV-l lr tT.*« pr>*i«- 

1. ir.itjtr r. . j.c ©j t«tti 

12 ;«ts •; loct*.«JcU . : ti*w : . th* t »st 

11 : c*n»?r ' i rr* ccnnrriny ir forr*ti-:s 

14 a <k> :©rp rrjr.?#r 

15 ;•«: Idv.tpolv.Attlv 
1< r#*l ifu; 

^ rft«* 

;i r • i; 

;4 •*. Irw^r* • ‘ ** fr. * ■ 
s j •: 

2 . • n**r k* kt:: re erruo*:*: xrF^.w:^ 

22 • a vaw.: :r r*** a failtf.i orrvmt 

2 i * a va u: cr t-<a* mrr km rJtimis ccrfOfitc kith Kt fa::.vf 2 

24 ds t-1 

35 c«ai:h;, » 0 

21 «r4 do 

2*> C- i*ll r. 

2 6 c«r**:s f • i • 1 
24 #r*4 do 

I# MS»t Vt KfK TC nrtF 7KI AAFA* -Vf I -S ; T^ZZ 
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n •. .** 'i * i 

32 • t.lt 

13 U-4 3 * ll.ii 

34 ? .»* ' 4 • • 2 2 31 

*5 t «<M 22 J* 

V * ;*# < • 3* -»'> 

V **v • . }» ); 

>4 • *,»# * * *f ; ‘ 

JJ nm» > • «' *4 

< : t ;*# ' ; ;. • t' 

4 . 1 !** . 2 • 11 

♦* * ' - d « 

*2 t :•« ;> * 2> 2 
44 24 - VC tl 

4 5 i.14 ;• * 15 V*. ' 

4 6 ;:•% ;< # 13'.>.*5 
4- r i»r 1- 1 132.31 
4* ?w 14 . yoo.ro 
4: ti*t il'.i ‘05.C0 

30 UM 2 ' - 20 . cc 

51 . th» call t© thi* rotroxitir.f :h 4 c<t fcr err'}-* ir. *>r «;u/ 

V 1 rch'scfc? r. :*uei 

t 2 op<n :c,ru*c« % t of :i« ut * j 

34 pfr »2 • 2.23 

3! writ***.*' «nt*i* win val j* 

54 r««d <!.♦* t»in 

5" write f«.« ‘ •r.t«4 v*.^« 

%4 rradi •: t*** 

44 V • < t*4X- / ICC 

4C t».n • tnur, • 4*. 

c: <»c 

€2 :».r. • t»;tt * 4*. 

4 3 c*l 1 plv^rj c«r.i <t.**4 ur.irr n.tr.n rot 3v . * ' M v 

<4 writ* 'll. 232 ftptlv . *cfclv 

45 1C3 UrMt 4# 1 6 . * , 

44 «:.J dr 

r 

4 4 t;.i 

WeibuU Plot With Confidence interval 

Profilc-ldu.-Iihoodbased confidence intervals can he determined for any percentile of interest using 
main program* similar to source code listings 5 and 6. The variable "pfrac" determines the percentile to 
be investigated. For example, using "pfrac ■ 0.30*' one can calculate data for the 30 percent life. Collec- 
tions of such confidence intervals were calculated, and the results arc provided in figures 7 and 8. These 
two figures are plotted using Wcibull coordinates so that the Wcibull cumulative distribution frequency 
will plot at a straight line. The test data points arc plotted at the positions of exact median ranks (Jaquelin 
1993 > 
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7.-*feibuti plot lor date of tatofa I with profNo-Wioljhootf-bMod conManca 
•ntorvolo. In*arv*a ahown art pofcttwfaa OO-parcant conManca Interval* of paroontila 

aaOmataa. 
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WgxtO. W ato u n ptof for data of tafcla I w«h proMHlMfihood baaed c onMan c a 
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Summary of Results 

Software ha* been developed to analyze fear faiifuc test data. The source code listing* discussed 
provided examples illustrating how to define the test data to be analyzed and how to call the appropriate 
subroutines to analyze the data. The software can be used to determine 

3 . Maximum likelihood estimates of the Wcibull distribution 

2. Data for contour plots of relative likelihood for two parameters 

3. Data for contour plots of joint confidence regions 

4. Data for the profile likelihood of the Weihull distribution parameters 
3. Data for the profile likelihood of any percentile of the distribution 

6. Profile-iikelihood'hased confidence intervals for parameters and/or percentiles of the distribution 

Contour plots of relative likelihood (not illustrated in this document) can be obtained from the output 
files as pros ided by source code listings 3 and 4. Profile-likelihood calculations for Weibull distribution 
parameters were not illustrated directly. A point for the profile-likelihood plot of the stupe parameter can 
be determined using subroutine h / ratio (appendix B). A profile likelihood for the scale parameter can be 
obtained using a source code in the manner of listings 5 and 6 and analyzing for the 63.2 pcreemiic 
fpfrac"» 0.632) 
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Appendix A 

Verification of Software 


The software was verified by reproducing results published by Meeker and Escobar H 998). The main 
program* to produce the figures in this appendix (figs. 9 to 14) were similar to the source code listings 
found in the main body of this document. 



•* 10/) 10* 10,4 10 * 10* 


* 

Fipurft %j— Contour* of oqu* W h a W h ood frtoi WibU> 9h *i$ pwffHf $ • 1 /a 
and Wata* pamtUr^odxp Inthtmannf offigwt S.1 from Matkar 
•nd Eacobar 



•J 100 10.2 10.4 10* 10* 


* 

ftfur* 10— Comourt at JoM conftdraa Mi. Wt*tou« *apt paramaitr 

0 a 1 /a and Wdbi*ac^para m alari| «•*?(»). In dttmannarolOfu* M 

•rom Moafcor and Eacobar (1 OH). 
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^® uf * P»a9jt» <Na9Hood it wthod to datarmlna eonfldonca intarvaf ter teal* 

parameter. Wrtteu* aeate paramotar^ ■ a*p <p). In tha marmar of ligura 6.9 Item 
•"•afcar ana Eaeotoar (1990). 



figunlZ P refiM ■i aOhoodmalhodtodi«afimnaconfldanoalnfiarv*teri*>apa 
pa ra matar. W afc t* #tapa parem at a r 0 » Mm. tn tea manor ot figure 0.0 (tew 
M aa fcar and Eaootoar (1990). 
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0000 0000 10000 12000 14000 16000 10000 20000 

io*perc#m Hit, to.i 

figure 1 3^CoXOMfi of equal |oint oonfiden ce level* for WMtouO shop# parameter 
and 10-percent We, fo,v WefeuR «hape paramet e r £ * 1/a. In the manner of 
figure 6.7 from Meeke r and Eecober (1990). 



10-percent Ha, *o.i 

figure 1 C ^roiie KMkhood method to determine c onidance Interval for 10-percent 
Ida, * 0 , 1 . m the manner of flgure 9.0 from Meeker end Eecober (1996). 


!<> 
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Appendix B 

Subroutine Descriptions and User Instructions 

SabfMdM CCHECK 

A subroutine to do "error*' checking on censoring array, ail values should be either I orO integers. 

Usage 

subroutine re heck in, censor) 

Arguments 

n (integer) sire of array '•censor" (input) 

censor (integer array) array of dimension (nj. censoring information (input) 


If all element* in airay "ten tor' equal either 0 or 1 . no action is taken. 

If any element does not meet the criteria, a line of text i* written to FORTRAN units numbered 6. 9. 
and 10. and then execution of the main program is halted. 

Other mbroutJnr* called 

None 
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Mrcariw CCHECK2 


A subroutine to do error checkin* on censoring array. All valuer should be either I or 0 integer* 

I'Mfr 

subroutine cchcck2 (n. censer) 

Argument* 

n < integer) size of array "censor" (input ) 

censor (integer arravt stray of dimension In), censoring information (input) 

Comment* 

If all elements in array "cer.soi" equal either 0 or I . no action is taken. 

Il any element does not meet the criteria, a line of text is w ritten to FORTRAN units numbered 6 and 
10, and then execution of the main program is halted. 

Other mbrmitfnm called 

None 
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MKMttoc nXSHAPEMAXLL 


A subroutine to calculate the maximum likelihood estimates of the scale parameter for a two- 
perameter Weibuil distribution. for a presumed shape and censored data. 

Usage 

subroutine fnshapemaxfl (times, censor, n. shape, scale) 

Ar g ume nt * 

times 
n 

censor 
shape 
scale 

C w nnscHt* 

This routine is based on the code published by Keats and Lawrence (1997). Although arguments 
arc pasted as single precision, calculations within the subroutine are done in double precision. 

For array “censor," a value of 0 indicates a test to failure whereas a value of 1 indicates a test sus- 
pended without failure (type 1 censoring). 

Other subrenting fated 

None 


(real array) array of dimension |n). test times (input) 

( integer) sire of array "censor" i input i 

(integer array) array of dimension (n). censoring information (input) 

( real i presumed Weibuil shape parameter ( input ) 

(real ) calculated Weibuil Kale parameter (output) 
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Subro utine PLVFEJICENT 


A subroutine to calculate a point on a profile-likelihood curve for a Wcibull distribution at a particu- 
lar cumulative distribution function (CDF) percentile of interest and a presumed value for the time to 
failure at that percentile. This routine can be called repeatedly using different values for the assumed 
time to failure tatimei to generate data for a profile-likelihood curve. 

I'mc* 

subroutine pitperrenr (time, censor, n. atime. pfrac. ratio, spblv. scblv) 


Argum e nt* 


lime 

(teal array’) 

array of dimension (n], test times (input) 

censor 

(integer array) 

array of dimension |n). censoring information (input) 

n 

(integer) 

size of array "censor" (input) 

ittmc 

(real) 

assumed time to failure (input) 

pfrac 

(real) 

assumed percentile of CDF corresponding to 



value of atime (input) 

ratio 

(real) 

profile ratio for best likelihood value [blvj (output) 

ipblv 

(real) 

shape factor corresponding to blv (output) 

icWv 

(real) 

wale factor corresponding to Mv (output) 

Comm^mU 




For array "censor.*' a value of 0 indicates a test to failure whereas a value of 1 indicates a test sus- 
pended without failure (type 1 censoring). 

Data are rescaled within the subroutine to avoid eweme values. 

The Weibull shape parameter resulting in the profile-likelihood value » found within the subroutine 
by an iterative method. The shape value is found to within a tolerance of 0.01 . A range for the shape must 
be provided to sun the iterative process, and the shape parameter resulting in the profile-likelihood value 
must be contained within that range. The subroutine uses a range from 005 to 8.0. These values could be 
adjusted to fit the needs of u particular data set. but one might encounter numerical problems for extreme 
Kale parameter values since it is used as an exponent in calculations. The routine writes a line of text to 
unit number 6 if the routine fails to converge within 1000 iterations. 


nmattt 

wiratKf 
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Mbrmrttec WBUNVERSE 


A Kibiouiine to calculate an inverse Weibul! cumulative dbtribuncm function (CDF) value. 

UMft 

Mibrouttne wMinvene < chape. Kale. frac. Iifci 

ArfMPrnti 


ahape 

•real) 

Wei bull shape parameter (input) 

Kale 

•real) 

Wcibull Kale parameter • input i 

free 

•real) 

Weihull CDF percentile of intcrot (input) 

life 

•real) 

time to failure for provided percentile of CDF 


CemmeaU 

None 


Oflwr wbrarttnn aMti I 
None 
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febcMCtue HTLSHAfE 


A subroutine to calculate Wcibull distribution likelihood-based profile data (likelihood ratk». 

V *age 

subroutine vilshapciz. censor, n. shapcmax. shapemint 

Arguments 

L 

censor 
n 

shapemax 
shapemin 

Comments 

For array "censorr a value of 0 indicates a test to failure whereas a value of I indicates a test 
suspended without failure (type I ccmoringi. 

No outputs are returned to the calling program. However* output is to FORTRAN unit 

number 9. if a file for unit number 9 was not yet opened, the user wj’j oe prompted for a filename. 

The likelihood profile ratio is calculated for 41 values of the shape factor* the values equally spaced 
between the values passed as shapemax and shapernm. The output written to unit number 9 consists of 
two columns of numbers of fixed length, w here the numbers are the shape factor and the calculated likeli- 
hood ratio 

Data arc rescaled w ithin the routine to avoid extreme value 

Other Mbrsatkws called 

* maxll 

fixnhapemaxll 

ratio 


(real array i array of dimension [n). test times (input) 

(integer anray > array of dimension [n], censoring information (input) 

(integer) size of array "censor** (input) 

(real) maximum shape value of range to he calculated (input) 

(real) minimum shape value of range to be calculated (input) 


NASA/TM-2O02 2UI0# 


26 


Subroutine WJLONE 


A subroutine to calculate Weibull distribution likelihood-ha ved joint confidence ratios and corre- 
sponding confidence levels to be plotted as contour plots. This routine presumes that the Weibull distribu- 
tion is parameterized by the shape parameter and 10-pcrccnt life. 

Usage 


subroutine n jloneii, censor, n. shapemax. shapemin. tenmax. tenmin) 


Arguments 

z 

(real array) 

censor 

(integer array) 

n 

(integer) 

shapemax 

(real) 

shapemin 

(real! 

tenmax 

(reab 

tenmin 

(real • 


array of dimension |n). test timer (input) 

aiTay of dimension (nj, censoring information (input) 

size of array ’’censor" (input) 

maximum shape value of range to be calculated (input) 

minimum shape value of range to be calculated (input) 

maximum 10-percent life to be calculated (input) 

minimum 10-percent life to be calculated (input) 


For array ~ censor," a value of 0 indicates a test to failure whereas a value of 1 indicates a test sus- 
pended without failure (type I censoring). 

Ko outputs are returned to the calling program However, output is written to FORTRAN units 
numbered 9 and 10. If files were not yet opened, the user will be prompted for filenames. 

The likelihood profile ratios and confidence levels an calculated for a full atny of 41x41 values 
(1681 lines) of the shape parameter and 10-percent life, the values equally spaced between the ranges 
defined by the values passed as shapemax to shapemin and ten max to tenmin, respectively. The output 
written to unit number 9 consists of three column* of numbers of fixed length, the numbers being the 
shape parameter, the 10-percem life, and the calculated likelihood ratio. The output written to unit number 
10 consuls of three columns of numbers of fixed length, where the numbers are the shape parameter, the 
10- percent life, and the calculated joint confidence level. 

Data are rescaled within the routine to avoid extreme values. 

NOTE: Confidence levels less than 10 percent are written to the output file as equal to 10 percent. 
Confidence levels greater than 99.99 percent are written to the output file as equal to 99.99 percent 

Otfwr m beau tines called 

w mat// 
n I ratio 

efudf (IMSL routine — see IMSI. ( 1997); IMSL it a registered trademark of Visual Numerics. Inc. » 


2i 
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Subroutine WJL TWO 


A subroutine to calculate WcibuII distribution likelihood-based joint confidence ratios and corresponding 
confidence levels to be plotted as contour plots. This routine pre turner that the Weibu’l distribution is 
parameterized by the shape parameter and scale parameter (or <ST2-perceni life). 

Usage 


subroutine u/Vrv.'o (z. censor, n, sbapemax, shapemin. scalemax. scale min) 


Arguments 

z 

(real array) 

Censor 

(integer array) 

n 

(integer) 

sbapemax 

(real) 

shapemin 

(real) 

scalemax 

(teal) 

seal em in 

(real) 


array of dimension [n]. test times (input) 

array of dimension |n), censoring information (input) 

size of array “censor" (input) 

maximum shape value of range to be calculated (input; 

minimum shape value of range to be calculated (input) 

maximum scale value of range to be calculated < input) 

minimum Kale value of range to be calculated (input) 


For array “censor,'' a value of 0 indicates a test to failure whereas a value of 1 indicates a test sus- 
pended without failure (type I censoring). 

No outputs an returned to the calling program. However, output is written to FORTRAN uniu 
numbered 9 and 10. If files were not yet opened, the user will be prompted for filenames. 

The likelihood profile ratios and confidence levels arc calculated for a full array of 41x41 values 
<1681 lines) of the shape parameter and scale parameter, the values equally spaced between the ranges 
defined by the values passed as sbapemax to shapemin and scalemax to scaJcmm, respectively . The out- 
put written to unit number 9 consists of three columns of numbers of fixed length, the numbers being the 
shape parameter, the scale parameter, and the calculated likelihood ratio. The output written to unit num- 
ber 10 consists of three columns of numbers of fixed length, where the numbers are the shape parameter, 
the scale parameter, and the calculated joint confidence level. 

Data are rescaled within the routine to avoid extreme values. 

NOTE: Confidence levels less than 10 percent are written to the output file as equal to 10 percent. 
Confidence levels greater than 99.99 percent arc written to the osaput file as equal to 99.99 percent. 


nmoxil 
wl ratio 

thtdf (IMSL routine— see IMSL (1997): I.MSL is a registered trademark of Visual Numeric*. Inc.) 
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A subroutine lo calculate the log likelihood value of a single life test Can be called repeatedly, with 
the outputs summed, to calculate the joint log- likelihood value of a data set 

L'«g* 

subroutine tlnghd c. sh. t, fail. Ilv) 

Arg u men ts 


ftC 

(renl> 

Weihul! scale parameter (input j 

fth 

(real) 

Weibull shape parameter finpuu 

t 

(real) 

test time (input) 

fail 

(integer ) 

flag indicating test to failure of test censored (input) 

Uv 

(real) 

log* likelihood value (output* 


Com menu 


For variable fail, a value of 0 indicates a test to failure whereas a value of I indicates a test suspended 
without failure (type I censoring i. 

Other subrout in es called 

None 
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A subroutine to calculate the numerator of a likelihood ratio of a data set unrig presumed Weibull 
shape and Kale parameters 

Uaafr 

subroutine wlratio aime. censor, n. shape, scale. !v) 


Arguments 


time 

censor 

n 

shape 

Kak 

Iv 


(real array > 

< integer array ) 
(integer) 

(real) 

(real) 

(real) 


array of dimension |n). test times (input) 

array of dimension jnj, censoring information (input) 

size of array "censor" (input) 

presumed Weibull shape parameter (input) 

presumed Weibull scale parameter (input) 

likelihood value, or value of numerator (output! 


For anas "censor." a value of 0 indicates a test to failure whereas a value of 1 indicates a test sus- 
pended without failure (type I censoring i. 

Other rt broe tlan catted 

H'foglt 
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wmaxll 


A subroutine to cakuUc (he maximum likelihood estimates of the teak parameter for a two- 
parameter Weibull distribution. for a presumed shape and censored data. 

Usage 


subroutine wmaxll fnmet. censor, n. shape, scale) 


Arguments 


times 

n 

censor 

shape 

Kale 


(real array) 
(integer) 
(integer array) 
(real) 

(real) 


an-ay of dimension (n), test limes (input) 

size of array ‘'censor'* (input) 

array of dimension [n). censoring information (input) 

calculated Wcibull shape parameter (output) 

calculated Wcibull scale parameter (output) 


This routine is based on the code published by Keats. Lawrence, and Wang (1997). Although argu- 
ments are passed as single precision, calculations within the subroutine are done in double precision. 

For array "censor,” a value of 0 indicates a test to failure whereas a value of 1 indicates a lest sus- 
pended without failure (type 1 censoring). 

The routine finds the maximum likelihood value for the shape parameter iteratively, where the shape 
parameter it the root of an equation. The routine is considered a* converged if the step size of the 
Newion-Raphson method is less than 0.0003. If the routine fails to converge after 100 iterations, o line is 
written to FORTRAN unit number 15. 

Other subroutines called 

None 
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